stats: silhouette score, ch, dbi

Part 1: stats

install.packages("wordspace") # contains dist.matrix function

also installing the dependencies ‘sparsesvd’, ‘iotools’

trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/sparsesvd_0.2.tgz'
Content type 'application/x-gzip' length 74751 bytes (72 KB)
==================================================
downloaded 72 KB

trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/iotools_0.3-2.tgz'
Content type 'application/x-gzip' length 154104 bytes (150 KB)
==================================================
downloaded 150 KB

trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/wordspace_0.2-6.tgz'
Content type 'application/x-gzip' length 2509428 bytes (2.4 MB)
==================================================
downloaded 2.4 MB

The downloaded binary packages are in
    /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//RtmpHMlGUN/downloaded_packages

library("wordspace")

Loading required package: Matrix

Attaching package: ‘Matrix’

The following object is masked from ‘package:flowCore’:

    %&%

dm <- dist.matrix(m, as.dist = TRUE)

Error: vector memory exhausted (limit reached?)

Saved the flowset object I created and load the object here - but I didn’t save the filtered df converted to matrix, which is needed for the stats below

fs <- readRDS("/Users/rhalenathomas/Documents/Data/FlowCytometry/PhenoID/Analysis/9MBO/prepro_outsjan20-9000cells/FlowSOMin.rds")

# note to Shuming - we should be able to go directly from flowset in preprocessing into this flowset but then the 9 samples (or however many will be in separate frames)
# I don't know where the actual values are located

# skip
#based on flowsom, the optimal k number of clusters = 5 for the sample (n=3000) 
metacl <- MetaClustering(fs$map$codes,
"metaClustering_consensus",
max = 15)

unique(metacl[fs$map$mapping[, 1]])

  dm = dist.matrix(m, as.dist = TRUE)

Error: vector memory exhausted (limit reached?)

# run the function to get the stats



krange <- 3:15 #range of number of clusters
li = stats(krange)  # doesn't contain silhouette 
# list of 2 

#silhouette score: ranges from -1  to 1 
#-1: bad clusters  0: neutral, indifferent  1: good clusters
#plot(krange, type='b', li[[1]][krange], xlab='Number of clusters', ylab='Average Silhouette Scores', frame=TRUE)

#Calinski-Harabasz index: 
# the highest value is the optimal number of clusters
plot(krange, type='b', li[[1]][krange], xlab='Number of clusters', ylab='Calinski-Harabasz index', frame=TRUE)

#Davies–Bouldin index: minimum score is zero
#the lowest value is the optimal number of clusters
plot(krange, type='b', li[[2]][krange], xlab='Number of clusters', ylab='Davies–Bouldin index', frame=TRUE)

The Calinski-Harabasz index show that K=8 is best number of clusters. But the Davies–Bouldin index shows lower is better. It is okay up until 9.

Try statistics on other clustering.

ch.pheno.k.50 = calinhara(m,df$phenograph,length((unique(df$phenograph))))  

Error in if (cln[i] < 2) cclx <- 0 : argument is of length zero

David Bouldin index for Phenograph

#calculate Davies–Bouldin index
# indiex.list, intracls, intercls 
# I'm not sure how this is working but it does run
dbi.phen.k.271 <- clv.Davies.Bouldin(cls.scatt.data(
      m,
      df$phenograph_clusterk271),
      intracls = "average",
      intercls = "average"
    )
  
dbi.phen.k.50 <- clv.Davies.Bouldin(cls.scatt.data(
      m,
      df$phenograph_cluster),
      intracls = "average",
      intercls = "average"
    )

make this calculations on the Seurat clustering

# read in the seurat object with the clustering values
seu <- readRDS("/Users/rhalenathomas/Documents/Data/FlowCytometry/PhenoID/Analysis/9MBO/prepro_outsjan20-9000cells/SeuratfromFlowsomPheno.rds" )

CH index and DBI from Seurat (Louvain) try one resolution first

dbi.louv.res.1 <- clv.Davies.Bouldin(cls.scatt.data(
      m,
      seu@meta.data$RNA_snn_res.1),
      intracls = "average",
      intercls = "average"
    )

Error in cls.id.vect.validity(clust, "clust") : 
  Bad usage: input 'clust' should be vector type.

NA

Calculate all CHI for Louvain

  ch_louv <-list()
l.range <- c("RNA_snn_res.0.1","RNA_snn_res.0.25","RNA_snn_res.0.5","RNA_snn_res.0.75","RNA_snn_res.1")
resolution = c(0.1,0.25,0.5,0.75,1)
prefix = "RNA_snn_res." 
# loop doesn't work because I can't index the seurat resolutions

ch = calinhara(m,seu@meta.data$RNA_snn_res.0.1,length((unique(seu@meta.data$RNA_snn_res.0.1))))
ch_louv["RNA_snn_res.0.1"] <- ch

ch = calinhara(m,seu@meta.data$RNA_snn_res.0.25,length((unique(seu@meta.data$RNA_snn_res.0.25))))
ch_louv["RNA_snn_res.0.25"] <- ch

ch = calinhara(m,seu@meta.data$RNA_snn_res.0.5,length((unique(seu@meta.data$RNA_snn_res.0.5))))
ch_louv["RNA_snn_res.0.5"] <- ch

ch = calinhara(m,seu@meta.data$RNA_snn_res.0.75,length((unique(seu@meta.data$RNA_snn_res.0.75))))
ch_louv["RNA_snn_res.0.75"] <- ch

ch = calinhara(m,seu@meta.data$RNA_snn_res.1,length((unique(seu@meta.data$RNA_snn_res.1))))
ch_louv["RNA_snn_res.1"] <- ch

Plot the CHI from Louvain clusters

Make a table with the CHI outputs

df.chi <- as.data.frame(Method, CHI.value)

Warning in as.data.frame.vector(x, ..., nm = nm) :
  'row.names' is not a character vector of length 7 -- omitting it. Will be an error!

Make another plot

#library(ggplot2)
ggplot(df.chi, aes(x= Method, y= CHI.value)) + geom_point() +
  theme(axis.text.x = element_text(angle = 90))

Part 2: plot

#transpose the csv so that seurat object has the right column and row
# the flow intensity values will be input as RNA expression
tm <- t(df2)
rownames(tm) <- colnames(df2)
colnames(tm) <- rownames(df2)


k <- 10 #change number of cluster here 
k <- 8 # best number of clusters according to CH index

metaClustering <- (metaClustering_consensus(fs$map$codes,k = k,seed=42)) # flowSOM clustering?

seu <- CreateSeuratObject(tm) # create a seurat object 
# check the seurat object by plotting some features
print(colnames(df2))
allAB <- colnames(df2)
VlnPlot(seu,features = c("CD56","AQP4")) # eerror in plot
DotPlot(seu, features = c("CD56","AQP4","CD24","GLAST","CD140a","CD29","CD184","CD71","O4","HepaCAM","CD133"))
# scale the data 
seu <- ScaleData(seu)

DotPlot(seu, features = c("CD56","AQP4","CD24","GLAST","CD140a","CD29","CD184","CD71","O4","HepaCAM","CD133"))
DotPlot(seu, features = allAB)
#do i need to normalize?
# seu <- NormalizeData(seu, normalization.method = "LogNormalize", scale.factor = 10000)

# create a UMAP
# to do so PCA must be run first 
seu <- FindVariableFeatures(seu) # needed for PCA or select all features
#dimentionality reduction methods: pca and umap

seu <-RunPCA(seu,seed.use = 42) # will use variable features by default
seu <- RunPCA(seu, features = allAB)

print(seu[["pca"]], dims = 1:5, nfeatures = 5) # tells you what AB are contributing to the 
seu <- RunUMAP(seu,features = allAB)

# Add the clustering data to the seurat object
seu <- AddMetaData(object=seu, metadata=metaClustering[fs$map$mapping[,1]], col.name = 'flowSOM.k.8')

DimPlot(seu, group.by = "flowSOM", reduction = "pca") 
DimPlot(seu, group.by = "flowSOM", reduction = "umap") # from normal PCA used as umap input it is mostly on blob, I get the same shape when using all AB features as PCA input

DotPlot (seu, features = allAB, group.by = "flowSOM")


#allow for color labeling
Idents(seu) <- "flowSOM"

FLowSOM clustering visually appear to be terrible.

I’ll just try the seurat clustering

# cluster using Louvain clustering
# using the square root of the number of cells for K might work better
#sqrt(73578) = 271.25
seu <- FindNeighbors(seu, dims = 1:10, k = 271)
seu <- FindClusters(seu, resolution = c(0.1,0.25,0.5,0.75,1))
clustree(seu, prefix = "RNA_snn_res.") + theme(legend.position = "bottom")
# this is  very slow - do not run again in notebook

library(clustree)
clustree(seu, prefix = "RNA_snn_res.") + theme(legend.position = "bottom")

# from Clustree 0.5 looks like the best resolution

# c(0.1,0.25,0.5,0.75,1)
resolutions = c("RNA_snn_res.0.1","RNA_snn_res.0.25","RNA_snn_res.0.5","RNA_snn_res.0.75","RNA_snn_res.1")

for (res in resolutions){
 print(DimPlot(seu, reduction = "umap", repel = TRUE, label = TRUE, group.by = res))
 
}

NA
NA

The clusters are not well separated, but better than with flowSOM (visually) Lets see how a dotplot looks

DotPlot(seu, group.by = "RNA_snn_res.0.5", features = allAB, scale = TRUE) # why is the intensity ploted?
# the data must not be working correctly


# maybe the object is not really correct at all
FeatureScatter(seu, feature1 = "CD56",feature2 = "CD24")
FeatureScatter(seu, feature1 = "CD29",feature2 = "CD184")

# maybe the data needs to be normalized
seu <- NormalizeData(seu)

DoHeatmap(seu, features = allAB, group.by = "RNA_snn_res.0.5")
# the groups look okay by heatmap. I'll try the different resolutions 

# c(0.1,0.25,0.5,0.75,1)
resolutions = c("RNA_snn_res.0.1","RNA_snn_res.0.25","RNA_snn_res.0.5","RNA_snn_res.0.75","RNA_snn_res.1")

for (res in resolutions){
 print(DoHeatmap(seu, features = allAB, group.by = res))
 
}

# try the dotplot again
# c(0.1,0.25,0.5,0.75,1)
resolutions = c("RNA_snn_res.0.1","RNA_snn_res.0.25","RNA_snn_res.0.5","RNA_snn_res.0.75","RNA_snn_res.1")

for (res in resolutions){
 print(DotPlot(seu, features = allAB, group.by = res, cols = c("blue","red")))
 
}

# what is happening?

Try feature maps

for (AB in allAB){
 print(FeaturePlot(seu, features = AB),slot = 'scale.data',min.cutoff = 'q10', max.cutoff ='90')
 
}

# autothreshold isn't great - I set quartiles 

#FeaturePlot(seu, features = "O4", min.cutoff = 1, max.cutoff = 25)
#FeaturePlot(seu, features = "O4", slot = 'scale.data',min.cutoff = 0.1, max.cutoff = 25)

#FeaturePlot(seu, features = "O4", slot = 'scale.data',min.cutoff = 'q5', max.cutoff ='99')

Try ridge plots

RidgePlot(seu, features = c("CD56","CD24"), ncol = 2) # missing values???
VlnPlot(seu, features = c("CD56","CD24"), ncol = 2) # also missing values???

I’ll save the seurat object and see about trying to cluster with phenograph. I believe I have read in the unalligned not normalized data

saveRDS(seu, "/Users/rhalenathomas/Documents/Data/FlowCytometry/PhenoID/Analysis/9MBO/prepro_outsjan20-9000cells/SeuratfromFlowsom.rds")

Phenograph clustering uses Jaccard coefficient between neruest neigbors sets and then id communities by louvain

# install
if(!require(devtools)){
  install.packages("devtools") # If not already installed
}
devtools::install_github("JinmiaoChenLab/Rphenograph")

library(Rphenograph)